In this notebook, we'll develop a model of a pandemic as it spreads in a susceptible population, and explore the effectiveness of interventions such as vaccination.
According to 2022 census records, the city of Banner Elk in North Carolina has a population of approximately 1,000 people. Banner Elk will be the center of our attention as we simulate the spread of an infection in the city. We introduce well-known models of infectious disease, the SI model, the Kermack-McKendrick (SIR) model, SIRS model, SIRV model and use them to explain the progression of the disease over the course of 2,000 days.
The following models were produced with the Epidemics on Networks (EoN) package for python
The following packages are necessary for this notebook:
import EoN
import networkx as nx
from collections import defaultdict
import matplotlib.pyplot as plt
import numpy as np
from PIL import Image
from collections import Counter
from IPython.display import HTML, display, Markdown, Latex
This model consists of nodes (individual) and edges (links between nodes).
Each node has an indepedent stage susceptible and infected. Any susceptible node(person) if infected, will again become susceptible and may be infected again at some point. Infection can take place at any time in this model.
Initially, we consider a model with no recovery rate. The susceptible compartment will flow entirely into the infected compartment over time.
This model is the most basic of all models in this document.
We can set parameters of the model as follows:
N = 1000 # Set the population
rho = 0.01 # Intial infected (10)
kave = 5 # Average degree
gamma = 0 # Recovery rate
tau = 0.002 # Transmision rate
S0 = (1-rho)*N
I0 = rho*N
SI0 = (1-rho)*kave*rho*N
SS0 = (1-rho)*kave*(1-rho)*N
t, S, I = EoN.SIS_homogeneous_pairwise(S0, I0, SI0, SS0, kave, tau, gamma, tmax=2000)
plt.plot(t, S, label = 'S')
plt.plot(t, I, label = 'I')
plt.xlabel('$t$')
plt.ylabel('Population')
plt.legend()
plt.show()
The SIR model is similar to the SIS model, but in this model after susceptible become infected they become R(t) or recovered. Where they cannot become infected again.
The susceptible equation is as follows:
dS/dt = -b * s(t) * I(t)
The recovered equation is:
dr/dt = k i(t)
The infected equation is:
ds/dt + di/dt + dr/dt = 0
These equations are utilized to develop our models below, we some variation.
N=1000 #population
gamma = 0.0045 #recovery rate
tau = 0.007 #transmission rate
kave = 5 #people we know
rho = 0.01 #initial infected
phiS0 = 1-rho
def psi(x):
return (1-rho)* np.exp(-kave*(1-x))
def psiPrime(x):
return (1-rho)*kave*np.exp(-kave*(1-x))
t, S, I, R = EoN.EBCM(N, psi, psiPrime, tau, gamma, phiS0, tmax = 2000)
plt.plot(t, S, label = 'S')
plt.plot(t, I, label = 'I')
plt.plot(t, R, label = 'R')
plt.xlabel('$t$')
plt.ylabel('Population')
plt.legend()
plt.show()
As seen above we follow a typical SIR Model.
The SIRS Model is a model based on the SIR, however, the recovered will eventually become S(t) again. This can add variation in any simulator, as persons can become infected again, after an infection has occured in a population region.
N = 1000
kave = 5 #expected number of friends/family
G = nx.fast_gnp_random_graph(N, kave/(N-1))
H = nx.DiGraph()
H.add_edge('I', 'R', rate = 0.0025) #Recovery rate
H.add_edge('R', 'S', rate = 0.001) #Recovered -> Susceptible
J = nx.DiGraph()
J.add_edge(('I', 'S'), ('I', 'I'), rate = 0.0065) #IS->II
IC = defaultdict(lambda: 'S')
for node in range(10):
IC[node] = 'I'
return_statuses = ('S', 'I', 'R')
t, S, I, R = EoN.Gillespie_simple_contagion(G, H, J, IC, return_statuses, tmax = 2000)
plt.plot(t, S, label = 'S')
plt.plot(t, I, label = 'I')
plt.plot(t, R, label = 'R')
plt.ylabel('Population')
plt.xlabel('$t$')
plt.legend()
<matplotlib.legend.Legend at 0x1e316070>
G = nx.grid_2d_graph(32,32) #each node is (u,v) where 0<=u,v<=32
#we'll initially infect those near the middle
initial_infections = [(u,v) for (u,v) in G if 14<u<18 and 14<v<18]
H = nx.DiGraph() #the spontaneous transitions
H.add_edge('Sus', 'Vac', rate = 0.0003)
H.add_edge('Inf', 'Rec', rate = 0.0025)
J = nx.DiGraph() #the induced transitions
J.add_edge(('Inf', 'Sus'), ('Inf', 'Inf'), rate = 0.01)
#calculate R for the graphs
def calc_R(v):
check_stat = []
R = 0
for time, w, *h in v:
check_stat.append([w])
c = Counter(i[0] for i in check_stat)
for x, y in c.items():
if(y > 0 and len(c.items()) > 0):
R += y
return R / len(c.items())
IC = defaultdict(lambda:'Sus')
for node in initial_infections:
IC[node] = 'Inf'
return_statuses = ['Sus', 'Inf', 'Rec', 'Vac']
color_dict = {'Sus': '#001aff','Inf':'#ff2000', 'Rec':'#008000','Vac': '#fffb00'}
pos = {node:node for node in G}
tex = False
sim_kwargs = {'color_dict':color_dict, 'pos':pos, 'tex':tex}
sim = EoN.Gillespie_simple_contagion(G, H, J, IC, return_statuses, tmax=2000, return_full_data=True, sim_kwargs=sim_kwargs)
display(Markdown(str(calc_R(sim.transmissions()))))
times, D = sim.summary()
newD = {'Sus+Vac':D['Sus']+D['Vac'], 'Inf+Rec' : D['Inf'] + D['Rec']}
new_timeseries = (times, newD)
sim.add_timeseries(new_timeseries, label = 'Simulation', color_dict={'Sus+Vac':'#E69A00', 'Inf+Rec':'#CD9AB3'})
sim.display(time=1, node_size = 4, ts_plots=[['Inf'], ['Sus+Vac', 'Inf+Rec']])
ani=sim.animate(ts_plots=[['Inf'], ['Sus+Vac', 'Inf+Rec']], node_size = 4)
ani.save('SIRV_animate.gif', writer="pillow", dpi=100)
1.46875
G = nx.grid_2d_graph(32,32) #each node is (u,v) where 0<=u,v<=32
#we'll initially infect those near the middle
initial_infections = [(u,v) for (u,v) in G if 14<u<18 and 14<v<18]
H = nx.DiGraph() #the spontaneous transitions
H.add_edge('Sus', 'Vac', rate = 0.0005)
H.add_edge('Inf', 'Rec', rate = 0.0025)
J = nx.DiGraph() #the induced transitions
J.add_edge(('Inf', 'Sus'), ('Inf', 'Inf'), rate = 0.01)
#calculate R for the graphs
def calc_R(v):
check_stat = []
R = 0
for time, w, *h in v:
check_stat.append([w])
c = Counter(i[0] for i in check_stat)
for x, y in c.items():
if(y > 0 and len(c.items()) > 0):
R += y
return R / len(c.items())
persons = {}
IC = defaultdict(lambda:'Sus')
for node in initial_infections:
IC[node] = 'Inf'
return_statuses = ['Sus', 'Inf', 'Rec', 'Vac']
color_dict = {'Sus': '#001aff','Inf':'#ff2000', 'Rec':'#008000','Vac': '#fffb00'}
pos = {node:node for node in G}
tex = False
sim_kwargs = {'color_dict':color_dict, 'pos':pos, 'tex':tex}
sim = EoN.Gillespie_simple_contagion(G, H, J, IC, return_statuses, tmax=2000, return_full_data=True, sim_kwargs=sim_kwargs)
display(Markdown(str(calc_R(sim.transmissions()))))
times, D = sim.summary()
newD = {'Sus+Vac':D['Sus']+D['Vac'], 'Inf+Rec' : D['Inf'] + D['Rec']}
new_timeseries = (times, newD)
sim.add_timeseries(new_timeseries, label = 'Simulation', color_dict={'Sus+Vac':'#E69A00', 'Inf+Rec':'#CD9AB3'})
sim.display(time=1, node_size = 4, ts_plots=[['Inf'], ['Sus+Vac', 'Inf+Rec']])
ani2=sim.animate(ts_plots=[['Inf'], ['Sus+Vac', 'Inf+Rec']], node_size = 4)
ani2.save('SIRV_animate.gif', writer="pillow", dpi=100)
1.4351585014409223
HTML(ani.to_jshtml())